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Context. We study the convection-pulsation coupling that occurs in cold Cepheids close to the red edge of the classical instability 
strip. In these stars, the surface convective zone is supposed to stabilise the radial oscillations excited by the /c-mechanism. 
Aims. We study the influence of the convective motions onto the amplitude and the nonlinear saturation of acoustic modes excited by 
/c-mechanism. We are interested in determining the physical conditions needed to lead to a quenching of oscillations by convection. 
Methods. We compute two-dimensional nonlinear simulations (DNS) of the convection-pulsation coupling, in which the oscillations 
are sustained by a continuous physical process: the /c-mechanism. Thanks to both a frequential analysis and a projection of the physical 
fields onto an acoustic subspace, we study how the convective motions affect the unstable radial oscillations. 

Results. Depending on the initial physical conditions, two main behaviours are obtained: (/) either the unstable fundamental acoustic 
mode has a large amplitude, carries the bulk of the kinetic energy and shows a nonlinear saturation similar to the purely radiative 
case; (ii) or the convective motions affect significantly the mode amplitude that remains very weak. In this second case, convection is 
quenching the acoustic oscillations. We interpret these discrepancies in terms of the difference in density contrast: larger stratification 
leads to smaller convective plumes that do not affect much the purely radial modes, while large-scale vortices may quench the 
oscillations. 

Key words. Hydrodynamics - Convection - Instabilities - Stars: oscillations - Methods: numerical - Stars: Variables: Cepheids 



1. Introduction 

The cold Cepheids located near the red edge of the clas- 
sical instability strip have an important surface convective 
zone that affects their pulsation properties (e.g. the reviews of 
Gautschy & Saio 1996; Buchler 2009). The first calculations, 
done without any convection-pulsation modelling, predicted red 
edges much cooler than the observed ones. In fact, predicting the 
red edge location requires a non-adiabatic treatment of this cou- 
pling as it was already stated by Baker & Kippenhahn (1965). 

Several time-dependent convection (TDC) models have 
been then introduced to address this coupling (Unno 1967; 
Gough 1977; Stellingwerf 1982; KuhfuB 1986; Xiong 1989). 
More recently, these different models have been largely 
improved and have succeeded in reproducing the cor- 
rect location of the red edge, despite their disagreements 
with the physical origin of the mode stabilisation (e.g. 
Bonoetal. 1999; Wuchterl & Feuchtinger 1998; Yecko et al. 
1998; Grigahcene et al. 2005; Dupret et al. 2005). 

However, these models rely on many free parameters (e.g. 
the seven dimensionless a coefficients used by Yecko et al. 
1998) that are either fitted to the observations or hardly con- 
strained by theoretical values such that they give almost identical 
results for different parameters sets (see also Kollath et al. 2002). 
Despite their own limitations (weak density or pressure con- 
trasts, for instance), 2-D and 3-D direct numerical simulations 
(DNS) are a good way to investigate the convection-pulsation 
coupling as the essential nonlinearities are fully taken into ac- 
count. The purpose of this paper is to present the results of two- 
dimensional DNS of the convection-pulsation coupling in which 
the acoustic oscillations are sustained in a self-consistent way by 



the /c-mechanism that is responsible for the radial oscillations of 
Cepheids (Muthsam et al. 2010). 

In Gastine & Dintrans (2008a) and Gastine & Dintrans 
(2008b) (hereafter Papers I and II), we have modelled the k- 
mechanism in Cepheids by a simplified approach, that is, the 
propagation of acoustic waves in a layer of perfect gas in which 
the ionisation region is shaped by a hollow in the radiative con- 
ductivity profile. We recall that the instability of Cepheid vari- 
ables relies on the blocking of the emerging radiative flux near 
the opacity bumps that are due to the ionisation of light elements 
like H or He. As the radiative conductivity is proportional to 
the inverse of opacity, a conductivity hollow therefore mimics 
such an opacity bump. By performing the linear stability of this 
simplified model, we have precisely determined in Paper I the 
physical conditions that are required to obtain unstable radial 
modes excited by /c-mechanism: ( i) the conductivity hollow must 
be sufficiently deep; ( ii) this hollow must be also located in a pre- 
cise zone inside the star, called "the transition region" (Zhevakin 
1963; Cox 1980), that defines the limit where the acoustic mode 
period n is close to the local thermal timescale r t herm as 



*P = 
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where the (c v To)Am term denotes the amount of internal energy 
embedded above the ionisation region and L is the star luminos- 
ity (the ratio between the two defining the thermal time needed 
for radiating this internal energy towards the surface). 

By starting from the most linearly-unstable setups, we have 
performed in Paper II the corresponding nonlinear study by the 
means of direct numerical simulations. Thanks to a powerful 
method that involves several projections of the computed fields 
onto suitable subspaces (e.g. Bogdan et al. 1993), both the am- 
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plitude and energy of each acoustic mode oscillating in the sim- 
ulation have been extracted. Their temporal evolutions have then 
emphasised the strong nonlinear coupling existing between the 
(unstable) fundamental acoustic mode and the (stable) second 
overtone. This preferential coupling relies on a 2:1 resonance as 
the fundamental period is twice the second overtone one (e.g. 
Simon & Schmidt 1976; Klapp et al. 1985; Buchler & Kovacs 
1986). It participates to the nonlinear saturation of the system 
and is also responsible for the well-known "Hertzsprung's pro- 
gression" observed in the luminosity curves of bump Cepheids 
(see Paper II for more details). 

Following these first two papers devoted to the purely ra- 
diative case, we next address in this work the possible influ- 
ence of convection onto the acoustic oscillations excited by k- 
mechanism. We have first slightly modified the initial physical 
setup to get a convective zone that overlaps with the hollow in 
radiative conductivity. We have then computed the correspond- 
ing 2-D nonlinear simulations and shown that a coupling be- 
tween the convective motions and pulsations can indeed occur 
(§2). Thanks to a frequential analysis (§ 3), we have next em- 
phasised the qualitative discrepancies between two specific sim- 
ulations that strongly differ from the role played by convection: 
in one simulation, the fundamental acoustic mode excited by k- 
mechanism seems to be not influenced by the underlying con- 
vective motions while in the other simulation, convection almost 
cancels it. 

In order to measure in more details the acoustic oscillations 
in these simulations, we have next applied a projection method 
similar to the one used previously in the purely radiative case. 
That is, we have determined the contribution of acoustic oscil- 
lations and convective motions to the energy budget (§4). The 
obtained results show again two opposite behaviours: (i) either 
the amplitude of the linearly-unstable acoustic mode is large and 
the nonlinear saturation seems similar to what was observed in 
radiative simulations; ( ii) or the amplitude of pulsations remains 
very weak and strongly modulated, while the bulk of kinetic en- 
ergy lies in convective motions. In this second case, the convec- 
tive plumes are quenching the oscillations in a similar way to 
what is expected to occur close to the red edge of the Cepheid 
instability strip. 

The last part of this study focus on the physical conditions 
suitable for the quenching of the oscillations by convective mo- 
tions (§5). We investigate the role of both the convective flux 
and the density stratification onto the mode amplitude and we 
show that a key parameter for the mode stabilisation could be 
the density contrast across the whole layer. Indeed, this stratifi- 
cation has an impact on the typical size of the convective plumes, 
as a weak (strong) stratification leads to large (small) vortices. A 
clear correlation between the integral scale of convection and the 
kinetic energy contained in acoustic modes is found, that sug- 
gests a screening-like effect of the mode pattern by convective 
plumes. We finally conclude in § 6 by discussing some inter- 
esting outlooks of this work in the direction of time-dependent 
convection theories. 



2. The model 

2.1. The convective setup 

Our system represents a local zoom around an ionisation region. 
It is composed of a 2-D layer (L x x L z ) filled with a monatomic 
and perfect gas with y = c p /c v = 5/3 and R* = c p - c v 
the ideal gas constant (c p and c v being the specific heats). As 
we are dealing with local simulations, the vertical gravity g = 




g 



Fig. 1. Sketch of our convective model. The blue curve corre- 
sponds to the radiative conductivity profile, the dashed line is the 
Schwarzschild criterion given in Eq. (6), and the large red arrow 
represents the radiative flux fbot that is imposed at the bottom. 
The red rolls correspond to the convective motions that develop 
in the layer middle. 



-ge z and the kinematic viscosity v are assumed to be constant. 
Following Papers I and II, the ionisation region is represented 
by a temperature-dependent radiative conductivity profile that 
mimics an opacity bump: 



K(T) = Kv 
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where 7b um p is the position of the hollow in temperature and <x, 
e and {A denote its slope, width and relative amplitude, respec- 
tively. We assume both radiative and hydrostatic equilibria, that 
is 



dpo 

dz 
dT 

dz 



-pog, 



(4) 



where is the imposed bottom flux. Following paper I and 
II, we choose the vertical dimension L z as the length scale, i.e. 
[x] = L z , the top density p top and top temperature 7\ op as density 
and temperature scales, respectively. The velocity scale is then 
(c p T top ) l/2 while time is given in units of [t] = L z /(c p T top ) l/2 . 

As our aim is to investigate the convection-pulsation cou- 
pling, we must ensure that convection will develop around the 
conductivity hollow. It means that our model should satisfy to 
the well-known Schwarzschild's criterion for the convective in- 
stability given by 
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Using Eqs. (4), it defines the following minimum value for the 
imposed bottom flux Ft, t 



Fbot ^ 



gKo(T ) 



(6) 



This criterion thus implies either a decrease in the gravity g or 
an increase in the imposed bottom flux Fbot- However, we have 
shown in Paper I that one crucial criterion needed to obtain un- 
stable modes excited by /c-mechanism is that the ratio of a local 
thermal timescale over the mode period is approximately (9(1), 
this is the relation (1) already given in the Introduction. Because 
of this timescale criterion, the most favourable setups were found 
to require large values in the gravity g to obtain a sufficient strat- 
ification (see Paper I for further details). As a consequence, the 
only acceptable way to have a convective zone in our model - 
while keeping the efficiency of /c-mechanism- is not decreasing 
the gravity g but rather increasing the bottom flux such that to 
satisfy Eq. (6). 

Figure 1 is a sketch of our convective model. The blue 
curve denotes the radiative conductivity profile from Eq. (2) 
that fulfills Schwarzschild's condition (Eq. 6, represented by the 
dashed black line), leading to the appearance of a convective 
zone shaped by the large red rolls where Kq(Tq) < c p F^ oi lg. 
Nevertheless, the need to increase the bottom flux to get a con- 
vective zone leads to a numerical difficulty that was not present 
in the pure radiative case. Indeed, we can roughly estimate the 
typical velocity w con v of a convective element by using a mixing- 
length argument of the form 



^bot 
PO 



1/3 



(7) 



With the typical values used in the former papers without con- 
vection, that is p to p - 2.5 x 10" 3 and Fb t - 2 x 10" 2 , one thus 
gets u com ~ 2 and therefore a Mach number Ma = u/c s > 1, 
with c s = ^|yR*T and T top = 1. It means that the parameters 
used in these radiative setups cannot be generalized to convec- 
tive simulations as they imply supersonic flows and shocks. 

To avoid such intricate shocks in the simulation, one should 
reduce the Mach number of the flow by both increasing the 
sound speed and decreasing the typical convective velocity. 
Increasing c s is done by simply doubling the value of the top 
temperature, that is, T top = 2. As Schwarzschild's criterion (6) 
imposes a large value of the bottom flux F^t, we then decrease 
the convective velocity w con v by increasing the top density to 
Ptop = 10" 2 , then a factor 4 compared to the one in the radia- 
tive case. We have also doubled the vertical extent of the layer 
to L z = 2 in order to satisfy the timescale criterion (1) while 
ensuring that the convective zone is deep enough. Finally, con- 
cerning the horizontal extent L x of the layer, the two following 
constraints remain to be taken into account: 

1. If the aspect ratio L x /L z is too small, the convective motions 
are too much constrained. Indeed, it is well known that a 
minimum horizontal wavenumber k x is needed to trigger the 
convective instability (e.g. Gough et al. 1976). As a conse- 
quence, dealing with small aspect ratio L x /L z can lead to a 
convective pattern too close to the Rayleigh-Benard one, that 
is, some large periodic rolls along the horizontal direction. 
In order to get a small-scale and as realistic as possible con- 
vective pattern, one should therefore consider domains with 
large aspect ratios. 

2. Nevertheless, dealing with a large aspect ratio has an im- 
portant numerical cost. Indeed, with a sixth-order finite- 
difference code as the one we are using, the mesh Reynolds 



number Re = uSx/v should not empirically be larger than 
~ 5. It means that for a given viscosity v and spatial resolu- 
tion, the horizontal extent is anyway limited. 

Within these limits, the aspect ratio is kept constant in every sim- 
ulation discussed in this work, with L x /L z = 4 for L z = 2 (thus 
Lx = 8). 

Concerning the parameters of the conductivity profile, the 
central temperature 7b um p has been adjusted to satisfy the crite- 
rion given in Eq. (1), that is *P = 1 at the location of the hollow. 
As in Papers I and II, its relative amplitude Ji is kept constant 
with the same value, that is - 0.7. At last, the slope <x and the 
width e of this profile are then set up to get as similar as possible 
profiles in the different DNS computed in this study. 

Table 1 summarizes the properties of the main numerical 
simulations performed to study the convection-pulsation cou- 
pling. The two coloured DNS, the G8 and G8H8 ones, are high- 
lighted because we will mainly focus on them in the follow- 
ing. The penultimate column of this table contains the value of 
the frequency 6l>oo of the fundamental unstable radial mode ex- 
cited by K-mechanism, that lies between 3 and 4 for every DNS 
computed in this study. The last column gives the value of the 
Rayleigh number that quantifies the strength of the convective 
motions and it is given by 



Ra: 



qL 4 

o conv 

YXCp 



(8) 



where L com is the width of the convective zone, x - Ko /poc p 
the radiative diffusivity and s the entropy. Tab. 1 shows that 
Rayleigh' s numbers of the DNS presented in this paper are close 
to 10 5 . This is well above the common critical values Ra c of this 
number from which strong convection is known to develop, e.g. 
Ra c ~ 10 3 for polytropic stratifications (Gough et al. 1976). 

2.2. Dimensional values of the physical quantities 

The numerical simulations computed in this work are performed 
in dimensionless units and can thus accommodate to a range of 
physical models. We can nevertheless redimension our units to 
determine if they are close enough to the realistic values of a typ- 
ical Cepheid star. We assume that the top of our layer is located 
under the surface at a temperature of 12000 K, it thus leads to 
the temperature scale [T] = 6000 K. For the G8 simulation given 
in Tab. 1, we have approximately T e [2., 13.], meaning that 
our layer covers a temperature range from 12000 K to 78000 K, 
that well encompasses the second Helium ionisation zone. For 
a 5 M Cepheid star with parameters close to the ones of 8- 
Cephei (i.e. R ~ 40 R Q , L ~ 2000 L Q and T eff ^ 6000 K), such a 
temperature range represents approximately 10% of the star ra- 
dius, leading to a length scale [z] = 1.9 x 10 9 m 1 . The density 
that corresponds to the temperature of 12000 K is approximately 
4 x 10" 5 kg/m 3 , giving a density scale [p] = 4 x 10" 3 kg/m 3 
for the G8 simulation. In such a star, one has Pbot/Ptop - 170, 
while our model has a lightly smaller density contrast of about 



Pbot/Pi 



top 



140. 



The timescale that corresponds to the G8 simulation can then 
been computed with the scaling dj ^c p T top , leading to [t] = 1 .3x 
10 5 s. It means that a frequency 6l>oo = 3.85 for the fundamental 
acoustic unstable mode corresponds in dimensional units to a 
mode period of 2.5 days. This value is of the same order than the 
period of £-Cephei, that is 5.36 days. 



1 The values have been taken in a 5 M© stellar model provided by the 
Helas network on http : //www . astro . up . pt/helas/stars/ . 
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DNS Gravity g Flux Fb ot Conductivity K m2LX 7b ump Width e Slope cr Viscosity v Frequency ojqo Rayleigh 



G8 


8 


4.5 x 10" 2 


io- 2 


6 


1 


1 


2.5 x 10~ 4 


3.85 


IO 5 


G8V5 


8 


4.5 x 10" 2 


io- 2 


6 


1 


1 


5 x IO" 4 


3.85 


5x 10 4 


G8H95 


8 


4.5 x 10" 2 


9.5 x IO" 3 


6.5 


1 


1 


5 x IO" 4 


3.84 


8x 10 4 


G8H9 


8 


4.5 x 10" 2 


9 x IO" 3 


7 


1 


1 


5 x IO" 4 


3.83 


IO 5 


G8H8 


8 


4.5 x 10~ 2 


8 x IO -3 


7.5 


1 


1.1 


5 x IO -4 


3.80 


2 x IO 5 


G7 


7 


4 x 10" 2 


io- 2 


5.5 


1.5 


0.8 


2.5 x IO" 4 


3.62 


8x 10 4 


G6 


6 


4 x 10" 2 


io- 2 


6 


1.5 


0.8 


5 x IO" 4 


3.35 


8x 10 4 


G6F7 


6 


3.7 x 10" 2 


IO" 2 


5.7 


1.5 


0.8 


3.5 x IO" 4 


3.36 


9x 10 4 


G6F7V4 


6 


3.7 x 10" 2 


IO" 2 


5.7 


1.5 


0.8 


4 x IO" 4 


3.36 


8x 10 4 


G6F5 


6 


3.5 x 10" 2 


IO" 2 


5.5 


1.5 


0.8 


3 x IO" 4 


3.36 


9x 10 4 


G6F5V4 


6 


3.5 x 10" 2 


IO" 2 


5.5 


1.5 


0.8 


4 x IO" 4 


3.36 


7x 10 4 



Table 1. Parameters of the numerical simulations done in this work. The two colored bold-typed ones emphasised the two simula- 
tions discussed in the following. For all these simulations, we assume 7\ op = 2, p top = 10" 2 , = 0.7 and L x /L z = 4. 



With these scalings, we can also check if the gravity and 
fluxes of the DNS are consistent with their stellar counterparts. 
The gravity g = 8 in our units then reads g = 0.84 m/s 2 . 
This value is very close to an estimate of the surface gravity 
g = GM/R 2 , that leads to g = 0.8 m/s 2 for £-Cephei. Concerning 
the fluxes, F hot = 4.5 x 10" 2 then becomes 5.2 x 10 8 J/s/m 2 . A 
good estimate of the real flux is obtained with F = L/4nR 2 , lead- 
ing for £-Cephei to 7.7x1 7 J/s/m 2 . Our simple model thus over- 
estimates the heat flux but in a rather moderate way with respect 
to what is usually done in global DNS (e.g. Brandenburg et al. 
2000). Such an accuracy in the heat flux is in fact inherent in 
the study of the K-mechanism as the thermal timescale has to be 
close enough to the dynamical one to fulfill the criterion (1). 

We can roughly estimate the kinematic viscosity of £-Cephei 
with the following law v ^ 1.2xl0~ 17 r 5 / 2 /p m 2 /s (see Chapman 
1954), that gives v = 2.5 x 10" 3 m 2 /s in the layer we focus 
on. A value vdns = 2.5 x 10" 4 in the DNS becomes vdns = 
6x 10 9 m 2 /s. We thus overestimate the viscosity by a factor 10 12 . 
Such an enormous value is a limitation of every DNS that can- 
not account for the small scale of convection because of the lim- 
ited numerical resolution (e.g. Chan & Sofia 1986). As a conse- 
quence, the viscosity is overestimated and leads to significantly 
lower Rayleigh and Reynolds numbers than the real ones. 

Table 2 sums up the differences between the stellar param- 
eters of a Cepheid of 5 M with the dimensional units of our 
local DNS. We notice that the values of our parameters are self- 
consistent as we both recover temperature and density variations 
that are compatible with the real ones, while keeping oscillation 
period, flux and gravity close enough to their stellar counterparts. 

2.3. Numerical methods 

As in Papers I and II, the linear stability analysis of the initial 
setup (4) has been computed thanks to the LSB code (Linear 
Solver Builder, Valdettaro et al. 2007). This spectral solver, 
based on the iterative Arnoldi-Chebyshev algorithm (Arnoldi 
195 1), gives both the complex eigenvalues X - r + icj and corre- 
sponding eigenfunctions if/ of a generalized eigenproblem of the 
form 



Ai/f = ABi/f, 



(9) 



where the sparse matrices A and B result from the projection 
of the differential operators on the Gauss-Lobatto grid. We es- 
pecially focus on the acoustic modes that are found to be lin- 
early excited by the /c-mechanism, that is those of which the real 
part of the eigenvalue is positive, r > 0. Furthermore, this linear 



a) 



Proc 1 



Transp. 



Proc 2 



b) 



Proc 1 



Proc 2 



Fig. 2. Sketch of the parallelisation of the ADI scheme, a) the 
initial domain decomposition of the Pencil Code for a two- 
dimensional domain in the (x, z) plane. The coloured lines de- 
note the data owned by the first processor (blue lines) and by the 
second one (green lines), b) the domain decomposition after the 
transposition needed to use the ADI solver. 

stability analysis has been computed without the consideration 
of the convective flux perturbations. This approximation, called 
the frozen-in convection, means that only the radiative contribu- 
tion of the total flux perturbation is allowed to change during the 
mode pulsation. We will see further that this approximation is 
well- suited in our case as the computed eigenfunctions overlap 
quite well with the vertical profiles extracted from the nonlinear 
simulation done with convection. 

Concerning the DNS, we have used again the Pencil Code 2 . 
This non-conservative code is an high-order finite-difference 
code (sixth order in space and third order in time) that con- 
serves the physical quantities up to the discretization errors of 
the scheme. On multiprocessor clusters, it takes advantage of 
the MPI libraries (Message Passing Interface) that allow to dis- 
tribute the computing tasks on several processors having their 
own memory. This code is fully explicit, except for the radiative 
diffusion term that is solved implicitly thanks to an Alternate 
Direction Implicit scheme (hereafter ADI) of our own already 
discussed in Paper I. However, this new study has needed an im- 
provement of this implicit solver as the convective motions that 
develop in our simulations require its parallelisation because of: 

1. the weakness of the mode growth rates r ~ 10" 3 computed 
from the new setups that imposes to carry out the simulations 
until a much longer time than in the purely radiative case. 

2. the convection itself that requires now a larger aspect ratio 
ensuring an efficient development of all additional spatial 



2 See http://www.nordita.org/software/pencil-code/ and 
Brandenburg & Dobler (2002). 
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Physical quantity 


DNS dimensionless values 


Scaling 


Length L z 


2 


d 


Ttop 


2 


Ttop 


^bot 




1 top 


Ptop 


10" 2 


Ptop 


Pbot 


1.4 


Ptop 


Oscillation period 


1.63 


d/ y]c p Ttop 


Gravity g 


8 


c p Ttop 1 d 


Flux 


4.5 x 10" 2 


Ptop(c p Ttopf 


Viscosity v 


2.5 x 10" 4 


d '\j c p Ttop 



DNS dimensional values Stellar values 



12000 K 12000 K 

78000 K 78000 K 

4 x 10" 5 kg/m 3 4 x 10" 5 kg/m 3 

5.7 x 10" 3 kg/m 3 7 x 10" 3 kg/m 3 

2.52 days 5.36 days 

0.84 m/s 2 0.8 m/s 2 

5.2 x 10 8 J/s/m 2 7.7 x 10 7 J/s/m 2 

6 x 10 9 m 2 /s 2.5 x 10" 3 m 2 /s 



Table 2. Comparison between the main physical quantities given in both the dimensionless units of the DNS and the corresponding 
dimensional values, and their stellar counterparts for a 5 M Cepheid. 



scales. Larger resolutions are thus compulsory to address the 
convection-pulsation coupling. 

In its 2-D version, the Pencil Code has a domain decompo- 
sition along the vertical direction such that the first direction of 
our ADI solver (the horizontal one) well fits in this topology. 
On the contrary, the second direction of the solver (the verti- 
cal one) is more tricky to parallelise as it needs a transposition 
of the data based on huge communications between processors. 
Figure 2 displays a sketch explaining how the parallelisation of 
the ADI solver is done. Once the tridiagonal system is solved 
in the second direction, one transposes once again the data to 
keep the original domain decomposition back. With this pow- 
erful algorithm, the semi-implicit nonlinear simulations of tur- 
bulent convection can then run on massively parallel clusters 
with the large resolutions induced by high Rayleigh numbers 
(see Tab. 1). Typically, our 2-D simulations of the convection- 
pulsation coupling were performed using a mid resolution of 
about 512 x512(i.e.512 grid points in each direction). 

2.4. 2-D DNS of the K-mechanism with convection 

With the parallel ADI solver for the radiative diffusion imple- 
mented in the Pencil Code (see Fig. 2), we advance in time the 
following hydrodynamic equations: 



( Dlnp 
Dt 



Du 
~Dt 



■ divw, 



1 



Vp + g + 2v (V • S + S • V In p) , 



(10) 



DT 1 

= — div KVT - (y - 1)7 divw + 2pvS 2 , 

Dt pc v 

where p, u and T denote the density, velocity and temperature, 
respectively. The operator D/Dt = d/dt + u • V is the usual total 
derivative, while S is the (traceless) rate-of- strain tensor given 
by 

1 Iduj duj 2 \ 

Finally, we impose that all fields are periodic in the horizontal 
direction, while stress-free boundary conditions (i.e. u z = and 
du x /dz = 0) are assumed for the velocity in the vertical one. 
Concerning the temperature, a perfect conductor at the bottom 
(i.e. flux imposed) and a perfect insulator at the top (i.e. temper- 
ature imposed) are applied. 



In order to ensure that both the nonlinear saturation and 
thermal relaxation are well achieved, the simulations were per- 
formed until very long times, typically t > 3000. With the eigen- 
frequencies computed from the linear stability analysis, that is 
d>oo ^3-4 (see Tab. 1), this approximately corresponds to 1500 
oscillation periods of the fundamental unstable acoustic mode. 

Figure 3 displays a snapshot of the modulus of the vorticity 
field (i.e. |V X u\) at a given time for the two simulations em- 
phasised in Tab. 1, the G8 (Fig. 3a) and G8H8 ones (Fig. 3b). 
The vorticity field highlights the convective motions that are, as 
expected by the sketch in Fig. 1, superimposed to the radiative 
conductivity hollow located approximately in the middle of the 
layer. Moreover, one notes that the typical size of the convective 
plumes differs in the two simulations as the eddies have a smaller 
scale in the G8 simulation than in the G8H8 one. Accordingly, 
the overshooting of convective elements into the bottom stably 
stratified layer is also more important in the G8H8 case than in 
the G8 one. It suggests a lower value of the Peclet number in 
that case as a faster thermalisation of the convective plumes at 
the interface between the convective/radiative regions is associ- 
ated with a larger penetration extent (Dintrans 2009). 

Beyond the qualitative differences seen in the vorticity field, 
a good way to compare the DNS is to study the temporal evolu- 
tion of average quantities, such as the vertical mass flux (pw z ), 
where the brackets (• • • ) denote a global average. Indeed, as we 
are doing simulations where both convective motions and os- 
cillations of acoustic modes are present, it is instructive to use 
a simple diagnostic that roughly separates their relative con- 
tributions. As the convective plumes have both ascending and 
descending movements, the average vertical mass flux removes 
quite well their contribution and one then gets a good tracer of 
the amplitude of the acoustic modes only. Figure 4 shows the re- 
sulting temporal evolution of (pu z ) for the G8 simulation (blue 
line) and the G8H8 one (green line). As displayed in the zoom 
in the bottom left corner, the time evolution shows an oscillating 
behaviour in both cases due to the radial oscillations of the fun- 
damental (and unstable) acoustic mode excited by /c-mechanism. 
In the G8 simulation, the amplitude experiences first an expo- 
nential growth until reaching the nonlinear saturation regime for 
time t ~ 1000. One notes that the transient duration is well com- 
patible with the growth rate of this mode that is about r ~ 10" 3 . 
At first glance, the temporal evolution of (pu z ) in the G8 simula- 
tion is similar to what has been observed in the purely radiative 
simulations of Paper II, that is, a linear growth of the amplitude 
before a saturation on a finite value on a timescale oc 1 /r. 

On the contrary, the dynamics of the G8H8 simulation radi- 
cally differs as the amplitude remains low compared to the radia- 
tive case and is highly modulated in time. As a consequence, no 
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Fig. 4. Temporal evolution of the mean vertical mass flux (pu z ) 
for the two simulations G8 (solid blue line) and G8H8 (solid 
green line). The two vertical dashed black lines define the bound- 
aries of the zoom displayed in the bottom left corner. 



clear nonlinear saturation is observed in this case and this kind 
of behaviour has not been found in the purely radiative models. 
It means that the acoustic oscillations are more influenced by 
convective motions in the G8H8 simulation than in the G8 one. 
To study further the physics involved in this coupling, we next 
perform a frequential analysis of these simulations. 



3. Frequential analysis 

3.1. Fourier decomposition 

An efficient tool to determine the modes that are present in a nu- 
merical simulation consists in taking first a double Fourier trans- 
form in space and time of the vertical mass flux 



pu z (x, z, t) 



pu z (k x ,z, 0)), 



where k x = (2n/L x )£ denotes the horizontal wavenumber, with 
£ = [0, 1, 2, . . .], while oj is the frequency. Second, one plots 
the power spectrum of pu z in the plane (z, co) for each value of 
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Fig. 5. Left : temporal power spectrum of the vertical mass flux 
pu z in the (z, cS) plane for the G8 simulation. Right : the result- 
ing spectrum after integrating over depth. In the £ = plane, the 
dashed blue vertical lines mark the position of the frequencies of 
the overtones obtained with the linear stability analysis. 



£ and "shark fin profiles" emerge about definite frequencies cor- 
responding to eigenmodes (Dintrans & Brandenburg 2004). The 
last operation requires an integration over the vertical direction 
z to get the final mean spectra for each £ of the initial field pu z . 

Figure 5 displays the results of this Fourier analysis applied 
to the G8 simulation, with the power spectrum of pu z for the first 
three values of £ on the left side, and the corresponding mean 
spectra on the right one. Many discrete peaks about given fre- 
quencies are clearly noticed in the radial plane £ = but also 
in the non-radial one £ = 1 , even the amplitudes are weaker by 
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Fig. 6. Left : temporal power spectrum for the vertical mass flux 
pu z in the (z, cS) plane for the G8H8 simulation. Right : the 
resulting spectrum after integrating over depth. 



an order of magnitude in this last case. As expected, the funda- 
mental radial acoustic mode ojqo = 3.85, which is the only one 
to be unstable linearly, has the largest amplitude in the £ - 
plane. This signature is seconded by numerous weaker ampli- 
tude peaks, of which the frequency is a multiple of ^oo- In fact, 
these peaks do no correspond to overtones as they do not overlap 
with the theoretical eigenfrequencies sorted out by the linear sta- 
bility analysis (the computed overtones are displayed as vertical 
dashed blue lines in Fig. 5). They rather correspond to harmon- 
ics of the fundamental mode and this is a signature of the large 
amplitude of this mode. Indeed, the presence of high- amplitude 
harmonics implies that the mode has a sufficient amplitude for 
generating them through a nonlinear cascade. To a certain extent, 
the same behaviour is obtained for the non-radial planes £ - 1 
and £ = 2, but with a much lower amplitude. 

The temporal evolution of the mean vertical mass flux dis- 
played in Fig. 4, which is similar to the ones obtained in the 
purely radiative case, could indicate a nonlinear saturation of 
the G8 simulation based on a 2:1 resonance between an un- 
stable driving mode and a linearly stable overtone that acts as 
a slave mode (see paper II for further details). But the Fourier 
analysis in Fig. 5 emphasises that there is no overlap between 
the frequencies of the DNS and the overtones ones. It means 
that the coupling between different eigenmodes is not favoured 
in the G8 simulation and a resonance-like phenomenon is not 
the physical process responsible for its nonlinear saturation. 
Furthermore, the presence of numerous large- amplitude har- 
monics is a guess of a mono-mode saturation, such as the one ob- 
served for instance with Van der Pol's oscillators (e.g. Krogdahl 
1955; Takeuti & Aikawa 1981; Buchler & Goupil 1984). 

The same Fourier analysis applied to the simulation G8H8 
gives quite different results compared to the G8 ones (Fig. 6). 



One first recovers the signature of some acoustic mode in the 
(z, cS) plane, especially in the radial £ = one where the fun- 
damental mode is well visible around the frequency ojqq = 3.80. 
But the disappearance of the harmonics whatever the value of 
£ means that the corresponding amplitudes of eigenmodes are 
much lower than in the G8 case. And indeed, the amplitude of 
the fundamental mode in the Fourier space is weaker by an order 
of magnitude in the G8H8 simulation. 

The differences in amplitude obtained previously in the time 
series in Fig. 4 are thus confirmed by the results obtained in this 
Fourier analysis: a high amplitude of the unstable acoustic mode, 
observed in the G8 simulation, goes with numerous harmonics 
and significant frequential signatures, while a weaker amplitude 
of this mode, observed in the G8H8 simulation, is coupled with 
the disappearance of these harmonics. 

3.2. Harmonic analysis 

In order to visualize the structure of an eigenmode oscillating in 
a direct numerical simulation, hydrodynamicists have developed 
a diagnostic tool called the "harmonic analysis". This method 
has been for instance used to visualize internal wave attractors 
propagating in a non-separable container (e.g. Hazewinkel et al. 
2008; Grisouard et al. 2008). The idea is to filter out the vertical 
velocity u z at a given frequency oj following the relation 

Qo>(x,z) =ff 2 u z (x,z,t)e ioj{t - h) dt, (12) 

where h - h + nT, with n an integer, and T = 2n/oj is the 
period. The obtained 2-D (complex) field U thus gives the pat- 
tern that evolves in time in the simulation like cos cot or sin o>t. 
In Grisouard et al. (2008), the chosen frequency was imposed 
by a forcing term of the form cos ojf r C t such that internal waves 
propagating in their simulation were necessarily at this single 
frequency ojf OYC (or its harmonics 2ojf OYC , 3ojf OYC , . . . ). In our case, 
the filtering frequency is of course given by the acoustic modes 
excited by K-mechanism and the resulting field U^x, z) should 
correspond to the eigenvectors that are solutions of the linear 
eigenvalue problem (9). We apply this method to the two simu- 
lations G8 and G8H8 for a filtering frequency equal to the fun- 
damental mode one, then oj = 3.85 and oj = 3.80 for the G8 and 
G8H8 simulation, respectively (see Tab. 1). The resulting mod- 
ulus \Ucj(x, z)\ is displayed in the figures 7 (G8) and 8 (G8H8). 

For the G8 simulation, the two-dimensional field |£/ Woo | 
shows only variations in the vertical direction and has a quasi- 
periodic behaviour in the horizontal one (upper panel in Fig. 7). 
This well corresponds to the pattern that we are expecting for a 
radial acoustic mode. Moreover, the convective plumes appear 
merely as very faint hints under the form of several small wig- 
gles. It means that the vertical motions at the frequency ojqo are 
dominated by the acoustic radial mode. A definitive confirmation 
comes from the comparison between the theoretical eigenvector 
and the horizontal average of | | given by 



l^oofe)! 



- _L f Lx 

L x Jo 



\U Um (x, z)\dx. 



(13) 



The result is shown in the bottom panel in Fig. 7 (solid black 
line). This vertical profile is then compared to the eigenfunction 
u z computed in the linear stability analysis (dashed blue line): 
the two lines overlap almost perfectly, meaning that the motions 
at the frequency co^ in the DNS have a vertical structure that 
corresponds to the unstable fundamental mode. 
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Fig. 9. Comparison between normalised vertical velocity (u z ) 
profiles for the fundamental mode (n = 0, £ = 0) according 
to the harmonic analysis of the DNS (solid black line) and the 
adiabatic eigenvector (dashed blue line). 



Fig. 7. Harmonic analysis of the vertical velocity field u z for 
the G8 simulation around the eigenfrequency ojqo = 3.85. Top: 
amplitude of the modulus of the filtered velocity |£/^ 00 (x,z)|. 
Bottom: normalised vertical profile of the horizontal average of 
the filtered velocity | | (solid black line) and the correspond- 
ing velocity eigenvector computed thanks to the linear stability 
analysis (dashed blue line). 



Harmonic analysis u = 3.80 
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Fig. 8. Harmonic analysis of the vertical velocity field u z for the 
G8H8 simulation around the eigenfrequency 6l>oo = 3.80. Top: 
amplitude of the modulus of the filtered velocity |£/ Woo (x,z)|. 
Bottom: normalised vertical profile of the horizontal average of 
the filtered velocity | £/^ 00 1 (solid black line) and the correspond- 
ing eigenvector of the linear stability analysis (dashed blue line). 

The same harmonic analysis performed on the G8H8 sim- 
ulation at the frequency ojoq = 3.80 leads to different results 
(Fig. 8). The 2-D modulus | Ucj 00 (x, z)\ indeed differs significantly 
from the previous one as large horizontal variations are observed 
and the expected invariance of the radial acoustic mode in that 
direction is broken (upper panel in Fig. 8). It means that con- 
vection has now an important impact on the mode pattern at 
this filtering frequency. The convective plumes strongly affect 
the acoustic oscillations in this DNS and this is consistent with 
both the temporal modulation of the vertical mass flux ampli- 
tude observed in Fig. 4 and the lower Fourier amplitude of the 
fundamental mode in Fig. 6. Despite this stronger coupling with 



convection, we nevertheless recover the structure of the unstable 
acoustic mode after averaging | | along the horizontal direc- 
tion (lower panel in Fig. 8). The agreement between the eigen- 
vector and the vertical profile is rather remarkable but not re- 
ally surprising as the average in the horizontal direction exactly 
amounts in the smoothing of the convective perturbations. 

The frequential analysis developed in this section (Fourier 
decomposition in the plane (z, co) and the harmonic filtering) 
confirms the discrepancies that were observed in the temporal 
evolutions of the mean vertical mass flux (pu z ) (Fig. 4): the 
acoustic oscillations are not influenced by convection in the G8 
simulation, whereas the convective motions have a strong effect 
on the modes evolution in the G8H8 simulation. We will now in- 
vestigate in more details this effect by the means of projections 
onto an acoustic subspace that give the temporal evolution of the 
kinetic energy embedded in the pulsations. 

4. Energy of oscillations 

4.1. The projection method 

We have already used a projection method onto an acoustic sub- 
space to determine the energy contained in the acoustic modes 
of radiative simulations (Paper II). This acoustic subspace was 
built from both normal and adjoint eigenmodes that were so- 
lutions of the linear equations for the perturbations. The need 
to consider the adjoint problem, and not only the regular one, 
was imposed by the non-hermiticity of the oscillations opera- 
tor. The radiative diffusion was indeed so large at the surface of 
the equilibrium models that strong non-adiabatic effects made 
the regular eigenmodes non-orthogonal (see also Bogdan et al. 
1993). However, in the current 2-D simulations with convection, 
the setup has significantly changed and the radiative diffusivity 
X oc 1 Ip is smaller than in the non-convective case (the top den- 
sity has been multiplied by a factor 4 in the convection setup, see 
§ 2. 1). The corresponding non-adiabatic effects at the surface are 
then weaker and we use a simpler approach than the one in Paper 
II. This method, based on the adiabatic eigenvectors, has been 
validated in Dintrans & Brandenburg (2004) by measuring the 
wave field generated by the oscillations of an entropy bubble in 
an isothermal atmosphere. 

To use this approximation, we must first ensure that the ver- 
tical velocity profiles in our 2-D DNS are sufficiently close to the 
adiabatic eigenfunctions of the linear problem for the perturba- 
tions. This is done in Fig. 9 where we compare the vertical pro- 
file obtained in the previous harmonic analysis of the G8 Simula- 
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tion (the solid black line in the lower panel in Fig. 7 and then also 
in Fig. 9) to the adiabatic eigenf unction u z (the dashed blue line 
in Fig. 9). The agreement between these two profiles is rather 
comfortable, indicating that we can shape our acoustic subspace 
from the adiabatic eigenvectors with a good confidence. This 
approach has the main advantage to simplify the computation 
of the kinetic energy contents as Lynden-Bell & Ostriker (1967) 
have demonstrated that adiabatic eigenvectors are mutually or- 
thogonal for this particular scalar product 



<^l,^2> 



fcpodz. 



(14) 



where the symbol f denotes the Hermitian conjugate. The acous- 
tic subspace is then given by the following relation 



ut(z, = Yj c tn(t)fan(z) with ($ t , n , , \jf €n ) = d w d nn > . (15) 

This defines the projection of the Fourier transform u^(z, t) of 
the velocity onto the adiabatic eigenvector ifr£ n (z) of degree £ and 
radial order n. The complex and time-dependent coefficient C{ n (t) 
entering in this projection is computed from 



Ct n (t) = tynlM 



= f P0*ln 

Jo 



• u t dz. 



(16) 



As an example, when a given eigenmode (€, n) is present in the 
simulation, this projection method leads to a projection coeffi- 
cient that behaves as ci n (f) oc e ia){nt , with a)£ n the mode frequency. 

4.2. The kinetic energy content 

The Parseval equality written in the Fourier space reads 



[ p u 2 (x,y,z)dV = I 

JV Jk 



p \u k \ 2 d 3 k, 



(17) 



where ul denotes the Fourier transform in space of u . In our case, 
we only perform a Fourier transform of the velocity field in the 
horizontal direction, therefore 



F tot 
^kin 



Jo 



p u 2 (x,z)dxdz= I y ypo\&t(z)\ 2 dz. (18) 



with still k x = (2n/L x )L With Eqs. (15-18), this leads to 

+oo +oo 

Oo = £2>„(oi 2 . (19) 



£=0 n=0 



This last equation means that the kinetic energy content of a sin- 
gle acoustic mode is simply equal to the square of its amplitude 
coefficient. In the 2-D simulations with convection, we have 



^kin 



(0 



^kin W + ^kin W, 



(20) 



where and E™™ are the kinetic energy embedded in 

acoustic modes and convective motions, respectively. In our 
setup, the fundamental radial mode (£ = 0, n = 0) is the only one 
that is excited by /c-mechanism therefore the acoustic energy is 
contained in low-degree and low-order modes (accordingly with 
the power spectra in Fig. 5). We can thus restrict the analysis to 
these modes only in the summations entering in Eq. (19), leading 
to 




Fig. 10. Temporal evolution of the acoustic kinetic energy ra- 
tio R(t) for the two simulations G8 (solid blue line) and G8H8 
(dashed green line) according to Eq. (22). 



4inW = £I>"(oi 2 + C„ nv (o. 

{=0 n=0 



(21) 



The ratio of the kinetic energy contained in acoustic modes can 
then be defined as 



3 5 

2Z |c <» w|2 



R(t) 



£=0 n=0 



(22) 



Figure 10 displays the temporal evolution of this ratio R(t) 
for the two simulations G8 (solid blue line) and G8H8 (dashed 
green line). In the G8 simulation, we see that the acoustic energy 
linearly increases until its nonlinear saturation above time f sat > 
10 3 , in a similar way to what has been observed with the mean 
vertical mass flux in Fig. 4. This timescale f sat is still compatible 
with the linear growth rate of the fundamental mode, that is, r ~ 
10" 3 and 4at ~ 1/t. Once this saturation is reached, the energy 
ratio remains large (i.e. > 70%, the remaining 30% being in the 
convection) and this behaviour is similar to the one observed in 
purely radiative simulations (see Paper II for further details). In 
other words, the acoustic oscillations are not much affected by 
the convective motions in this simulation. 

On the contrary, the acoustic energy ratio remains very weak 
in the G8H8 simulation. Indeed, despite some transient increases 
during which non-trifling values R 10% are obtained, the av- 
erage ratio is R < 1-5% and convective motions are responsible 
for the bulk of the kinetic energy content. In this case, the radial 
oscillations excited by /c-mechanism are thus quenched by con- 
vective plumes. This situation is relevant to the red edge of the 
classical instability strip, where the unstable acoustic modes are 
supposed to be damped by the surface convective motions. 

4.3. Mean-field analysis 

Beyond the temporal evolution of the energy ratio displayed in 
Fig. 10, it is interesting to precisely locate the zones where ki- 
netic energy is due to acoustic modes and convective motions, 
respectively. Towards this goal, a simple mean-field analysis is a 



10 



T. Gastine and B. Dintrans: Convective quenching of stellar pulsations 



valuable tool to separate the mean-field component of the motion 
(the acoustic oscillations in our setup) from its fluctuating part 
(here the convective plumes). Such approach have been used for 
instance to check the mixing-length theory of convection thanks 
to velocity correlations (Chan & Sofia 1987). 

In our simulations, we are interested in the mean-field veloc- 
ity correlations that can account for the acoustic modes and the 
convective plumes. We first separate the vertical velocity field 
into a mean part and a fluctuating one 



u z = (u z ) + u' 



(23) 



where (• • • ) is an horizontal average. We then square this expres- 
sion and average along the horizontal direction to get 



(u 2 z ) = (u z ) 2 + (u' 2 ). 



(24) 



This equation is the mean square velocity that can be split be- 
tween an acoustic contribution and a convective one, i.e. (u 2 ) = 
acoustic + convection. To separate these different contributions, 
we recall that (u z ) is a good proxy to measure the acoustic os- 
cillations, as the vertical convective motions almost vanish after 
their averaging in the horizontal direction (ascending and de- 
scending plumes are comparable). The acoustic contribution to 
Eq. (24) thus reads 



Acoustic G8 
Convection G8 
Acoustic G8H8 
Convection G8H8 




Height 

Fig. 11. Vertical second-order velocity correlations for the sim- 
ulations G8 (blue lines) and G8H8 (green lines). The acoustic 
contribution (Eq. 25) is displayed as solid lines, while the con- 
vective one (Eq. 26) is plotted as dashed lines. 



<Kcousfe) = < u z) 



(25) 



With Eq. (24), we then extract the contribution of the convective 
motions to the mean square velocity field as 



su z com ( Z ) = «> = «> - (u z y 



(26) 



Fig. 11 displays the vertical profiles of Sul cous and du 2 om 
for the simulations G8 (blue lines) and G8H8 (green lines). In 
both simulations, one notes that the maximum of the convective 
contribution (dashed lines) is approximately located at the layer 
middle, corresponding to the radiative conductivity profile dis- 
played in Fig. 1, while the acoustic part (solid lines) reaches its 
maximum close to the surface because of the eigenvector shape. 

However, the acoustic contribution to the velocity correla- 
tion for the G8 simulation is larger than the convective one and 
a profile similar to the velocity eigenfunction is obtained (see 
Fig. 9 for instance). On the contrary, for the G8H8 simulation, 
this acoustic part Sul cous is much more weaker than the one due 
to convective motions. 

These trends can easily be linked to what was obtained pre- 
viously with the projection formalism (see Fig. 10). Indeed, the 
vertical profiles in Fig. 1 1 can be integrated to get a ratio between 
the acoustic oscillations and the (total) mean square velocity (u 2 ) 
as 



^lim - 



f 

Jo 



^acous * 



Jo 



(27) 



)dz 



This last equation, derived from a mean-field analysis, contains 
the same physical information than in Eq. (22) as the obtained ra- 
tio is very close to the one computed with the projection method 
in the saturated regime (Fig. 10): R\{ m ^75% for the G8 simula- 
tion and R\ im 3% for the G8H8 one. The second-order correla- 
tions Sul cous and 6u 2 om are thus also a good tool to separate and 
locate the different contributions to the flow velocity. 

The amount of energy contained in convective motions and 
in oscillations has been studied thanks to both the projection 



formalism already used in purely radiative simulations, and the 
second-order correlations of velocity. One recovers the main 
physics underlined in the frequential analysis, that is, in the G8 
simulation the mode amplitude is strong and the oscillations 
clearly dominate the motions in the flow, in opposition to the 
G8H8 simulation in which convection is dominant. 



5. The mode quenching 

We are now going to focus on the physical process that lead to 
the different behaviours obtained in our different DNS. In fact, 
as presented in Tab. 1, the main physical parameters of these 
simulations are very close: for instance, concerning the G8 and 
G8H8 simulation, there is mainly a 20% variation of the mean 
radiative conductivity K max . Some of the parameters of the con- 
ductivity hollow (e.g. Tbump and cr) also differ as they have been 
adjusted to satisfy the criterion given in Eq. (1) in each DNS. 
We can classify the simulations given in Tab. 1 in three different 
categories: 



- First, the simulations similar to the G8 one, in which the am- 
plitude of the unstable acoustic mode is important (R > 70%) 
and that doesn't seem much affected by convection. This 
case concerns the simulation G7 of Tab. 1 . 

- The second category, similar to G8H8, contains the DNS, 
in which the amplitude of the acoustic mode is very weak 
(R < 10%). In this case, the oscillations are quenched by 
convective plumes that dominate the flow. This case con- 
cerns also the G6 and G6F7 simulations of Tab. 1 . 

- At last, there are intermediate simulations, in which the 
amplitudes of oscillations and convection are comparable 
(R ^ 30 - 50%). These simulations also show important 
temporal modulation of the amplitude of oscillations: i.e. the 
standard deviation of the ratio of energy is larger than in pre- 
vious categories. It mainly concerns the simulations G8H9 
andG6F5 of Tab. 1. 
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Tentfe) = C p {pU z T') + C p {pU z ) 0. 



(32) 
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Fig. 12. Mean vertical profiles of radiative T ra d (solid blue line), 
turbulent enthalpy T CO rw (dashed green line), kinetic Tkin (dotted 
black line), modes Tmod (dot-dashed magenta line) and total Ttot 
(dotted red line) fluxes for the G8 (a) and G8H8 (b) simulations. 
All fluxes have been normalised with respect to the bottom value 

5.1. The role of the convective flux 

A first attempt in explaining these different physical behaviours 
is to compare the heat transport between the various DNS that 
are summarised in Tab. 1 . Especially, the study of the respective 
amount of heat carried out by convection and radiation may lead 
to a natural explanation to the mode quenching observed in the 
G8H8 simulation as, from an intuitive point of view, the higher 
the convective flux the stronger the quenching. 

The vertical heat transport is ensured by the radiative Trad, 
enthalpy Tent an d kinetic Tkm fluxes given by (e.g. Hurlburt et al. 
1984) 



'Tmdfc t) 



- (K(T)VT) , 
c p (pu z T) , 



^Tkinfe t)=~ ((U 2 X + u\)pu z ) , 



(28) 



where the brackets denote an horizontal average. The enthalpy 
flux Tent can be developed in mean and fluctuating components 
from the following decomposition for the total temperature T 



T(x,z, t) = (T)(z, t) + T'(x,z, t), 

= T m (z) + 6(z,t) + T'(x,z,t), 



(29) 



where Ths is the hydrostatic background temperature profile, 6 is 
the temperature eigenfunction of the unstable acoustic mode and 
T is the turbulent fluctuating temperature around the horizontal 
average (T). The enthalpy flux then reads 

Tentiz, t) = c p (pu z r) + c p (pu z ) (T m + 0) . (30) 

We now average in time this last expression over a multiple of 
the mode period to get 



Tent(z) = Tentiz, t) = C p (pU z T') + C p {pu z ) 6 + C p T m {pU z ). (31) 

Finally, as there is no average mass flux over an oscillation pe- 
riod (i.e. (pu z ) = 0), the enthalpy flux carried by both the con- 
vection and acoustic oscillations is given by 



The first term of the right hand side denotes the turbulent trans- 
port of enthalpy usually defined as the convective flux (e.g. 
Hurlburt et al. 1984), while the second one corresponds to en- 
thalpy transported by the acoustic modes. This second contribu- 
tion is usually negligible in simulations of compressible convec- 
tion because of the weakness of modes that are driven by tur- 
bulent fluctuations in pressure (Bogdan et al. 1993). In our case, 
the modes can efficiently be excited through the /c-mechanism 
and their contribution to the heat transport should be taken into 
account. In order to determine the net contribution of the acous- 
tic modes to the total flux, we thus separate the usual turbulent 
enthalpy part due to the convective plumes (hereafter denoted 
Tconv) from the enthalpy transport due to modes (hereafter Tmod)- 
By taking the average in time of Eqs. (28) and using Eq. (32), 
one finally gets 



T md (z) = -(K(T)VT), 



T'comiz) = c p {pu z T'), 

Tydn(z) = \{(M 2 X + Uj)pu z y 



(33) 



v Tmodiz) = C p {pU z )6. 

The vertical profiles of these fluxes (normalised to the imposed 
bottom value Fb t) are given in Fig. 12 for the G8 (upper panel) 
and G8H8 (lower panel) simulations. For both DNS, the bulk 
of the total flux is transported by the radiative flux, except in 
the convective zone where T com reaches 15% of the total flux 
for the G8 simulation and 40% for the G8H8 one. We also no- 
tice a larger overshooting in the second DNS, as the downdrafts 
penetrating in the radiative zone are stronger (see Fig. 3). As 
expected, this significant penetration is associated with both a 
negative convective flux and a small value of the kinetic flux 
(~ 5%). The convective fluxes in our DNS are nevertheless rel- 
atively weak compared to what is expected to occur in Cepheid 
stars (see e.g. Yecko et al. 1998). These limited values are a di- 
rect issue of the rather moderate Rayleigh number numerically 
affordable in such kind of DNS. 

Concerning Tmod, one notes that it is hardly as important as 
the convective flux in the G8 simulation, while it is almost van- 
ishing in the G8H8 one. In the G8 simulation where the ampli- 
tude of the unstable acoustic mode is the largest, about 10% of 
the heat flux is carried out by acoustic modes. This quantity is 
thus another good tracer of the significance of the amplitude of 
the acoustic modes compared to the convective motions, as Tmod 
becomes significant when the amplitude of the acoustic oscilla- 
tions is strong enough. 

The differences in the amount of heat transported by turbu- 
lent convection may be a first hint on the physical origin of the 
mode quenching observed in our DNS. To confirm this hypoth- 
esis, we have thus gathered the values of the maximum of the 
normalised convective flux for all the DNS in Tab. 1. Figure 13 
displays the amount of energy contained in acoustic modes with 
respect to these maximum values. For moderate to stronger con- 
vective fluxes, we recover a significant relation between the frac- 
tion of T C om and the efficiency of the /c-mechanism, that is, the 
energy ratios in acoustic modes become smaller with increasing 
convective fluxes. But this correlation disappears in the region 
of smaller convective fluxes with T CO rw £ 20%. Indeed, one ob- 
serves here that the efficiency of the /c-mechanism is independent 
of the value of the convective flux as the modes may be either 
quenched or excited for similar values of f conv . 
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G8 


0.42 


76% 


G8V5 


0.42 


77% 


G8H95 


0.44 


71% 


G8H9 


0.47 


42% 


G8H8 


0.53 


3% 


G7 


0.44 


70% 


G6 


0.55 


4% 


G6F7 


0.50 


15% 


G6F7V4 


0.50 


10% 


G6F5 


0.48 


41% 


G6F5V4 


0.48 


28% 



Table 3. Values of the density scale at the layer middle accord- 
ing to Eq. (37) and ratio of energy contained in acoustic modes 
according to Eq. (22). 



Fig. 13. Ratio of the kinetic energy in acoustic modes displayed 
as a function of the maximum of the normalised convective flux 
^onv for the different DNS presented in Table 1 . The blue dots 
correspond to the mean kinetic ratio obtained with the projec- 
tion method, while the vertical gray bars represent the standard 
deviation. 




Fig. 14. Mean vertical profiles of density p for the G8 simulation 
(solid blue line) and the G8H8 simulation (dashed green line). 



Fig. 13 thus emphasises that the physical origin of the mode 
quenching does not solely rely on the amount of the convective 
flux, as the acoustic oscillations are also damped in DNS where 
the convective transport is inefficient. This result indicates that 
other underlying physical processes are involved in this mode 
quenching. 

5.2. The role of the stratification 

The differences between the physical parameters of the DNS 
given in Table 1 seem at a first look negligible. However the 
variations of gravity g, flux F\> ot and radiative conductivity K mSLX 
lead to significantly modified equilibrium fields, especially con- 
cerning the density profile. The mean vertical profile of density 
is close to the hydrostatic equilibrium, that reads 



dlnp 

dz 



"hot 



K(T) R* 



(34) 



where the radiative conductivity K(T) is given by Eq. (2). 
Increasing J3 = Fb 0t /K mSLX , or decreasing the gravity g then leads 
to a weaker derivative of the density. Fig. 14 displays the mean 
vertical density profile for the simulations G8 (solid blue line) 
and G8H8 (dashed green line). It shows that a 20% change of 
the radiative conductivity leads to significantly different density 
values, especially at the bottom of the layer where one notes a 
factor of two between the simulations. 

To compare the differences in the density stratification be- 
tween the simulations of Tab. 1 , we may use an analytical ap- 
proximation of the density scale H p = -(dlnp/dz)~ l at the layer 
middle. Starting from Eq. (34), one can approximate the temper- 
ature profile by 



dT 
~dl 



^bot 

K(T) 



T{z)^-f3{z-L z ) + T x 



top? 



(35) 



where we suppose that the mean temperature profile is linear, 
that is the conductivity profile is neglected and assumed to be 
K(T) ^ K max . The hydrostatic equilibrium then reads 



dlnp 

dz 



1 



(36) 



Lz) + ^top 

We can then approximate the density scale H p at the layer middle 
by 



H, 



P middle 



PLJ2 + r top 
p-glR* 



(37) 



This approximation allows us to compare the different simula- 
tions of Tab. 1 . The results given in Tab. 3 emphasises a correla- 
tion between the density profile and the ratio of energy contained 
in acoustic modes, that is the weaker H p the larger Rn m . On the 
contrary acoustic modes are quenched by convective motions if 
the value of H p increases (simulations G8H8 or G6 for instance). 

We can also extract the real values of H p directly from the 
DNS field in order to check that the results obtained in Tab. 3 
are confirmed without the constant radiative conductivity profile 
approximation done in Eq. (37). Fig. 15 displays the ratio R(t) as 
a function of the mean density scale (H p ) for the different DNS. 

Blue dots corresponds to R(t) (temporal average) and vertical 
bars to the standard deviation of R(t). We recover the correla- 
tion between the value of the density scale and the significance 
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Fig. 15. Ratio of the kinetic energy in acoustic modes displayed 
as a function of the mean density scale (H p ) for the different 
DNS presented in Table 1 . The blue dots correspond to the mean 
kinetic ratio obtained with the projection method, while the ver- 
tical gray bars represent the standard deviation. 



Fig. 16. Normalised kinetic energy power spectrum in the middle 
of the convective zone for the G8 simulation (solid blue line) 
and the G8H8 simulation (dashed green line). The kinetic energy 
spectrum E(k) ~ k~ 3 has been superimposed as a solid black 
line. 



of pulsations compared to the convective motions. What's more, 
Fig. 15 emphasises the three different behaviours presented be- 
fore, as one gets the two extrema (high amplitude of modes or 
mode-quenching by convection) but also "mixed"-simulations 
with R ^ 30 - 50% and very large standard deviations mean- 
ing that an important modulation of pulsations occurs. 

5.3. Spectrum and integral scale 

Fig 15 seems to indicate that the density stratification plays an 
important role on the stabilisation of acoustic modes excited by 
/c-mechanism. According to the mixing length theory (Vitense 
1953; Bohm- Vitense 1958), more important values of H p means 
that the typical size of convective eddies increases. This assump- 
tion is besides displayed on Fig. 3, where the differences be- 
tween the size of the vortices in the G8 and G8H8 simulations 
are noticeable. To determine the role played by the different 
scales of the flow onto the pulsations, we now focus on the power 
spectra of the different DNS. 

The power spectrum is obtained with a Fourier transform of 
the kinetic energy along the horizontal direction : 



E(k,z) 



Jo 



p(u 2 x + u 2 z )e ikx dx. 



(38) 



The normalised power spectrum at the middle of the convective 
layer (z - 1.2) is then displayed on Fig. 16 for the two simu- 
lations G8 and G8H8. One notes that the solid blue line, corre- 
sponding to G8, reaches its maximum to an higher value of k x 
than the dashed green line (G8H8). The energy spectrum is thus 
shifted to smaller scale in the G8 simulation than in the G8H8, in 
which the kinetic energy is stored in larger scales. It corresponds 
to what is observed on the vorticity fields of Fig. 3, where the 
convective eddies are larger in the G8H8 simulation than in the 
G8 one. 

We can also mention that we recover in the kinetic energy 
spectra a scaling law E(k) ~ k~ 3 . This result corresponds to the 
standard dimensional analysis of the enstrophy cascade in two- 
dimensional turbulence (e.g. Lesieur 1997). 




0.32 0.34 0.36 

Integral scale 



Fig. 17. Ratio of the kinetic energy in acoustic modes displayed 
as a function of the integral scale (Eq. 39) for the different DNS 
presented in Table 1 . The blue dots correspond to the mean ki- 
netic ratio obtained with the projection method, while the verti- 
cal gray bars represent the standard deviation. 



From Eq. (38), we can determine the integral scale, that is the 
scale associated with the most energetic structures of the flow. 



f 

Jo 



E(k,z)k~ l dk 



Jo 



(39) 



E(k, z)dk 



Once again, we compute the integral scale at the middle of 
the convective zone for the different DNS of Tab. 1. Fig. 17 
displays the value of R(t) with respect to the integral scale £ mt . 
The vertical gray bars emphasize the standard deviation of R(t). 
Fig. 17 illustrates the same trend as in Fig. 15: the larger the in- 
tegral scale, the smaller the energy ratio. The integral scale is 
in fact another way to get the information obtained with the H p 
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length scale, but relies on the properties of the turbulence of the 
flow. 

The comparison of these length scales in our DNS is a first 
hint on the physics responsible for the mode quenching. In the 
light of Figs. 15-17, the density stratification seems to have an 
important impact on how the oscillations develop in the pres- 
ence of convection. With the increase of the size of the eddies 
(increase of £ mX ), the acoustic part in the energy budget is gradu- 
ally vanishing. 

We may explain this mode quenching by the spatial distribu- 
tion of the convective plumes: large eddies entail an important 
"screening effect" onto the acoustic radial mode. One can imag- 
ine that the amplitude of such mode is affected by the motions 
that are not coherent to its vertical structure: in this case large ed- 
dies with significant amount of energy distributed both in radial 
and non-radial motions is incoherent with the "main" velocity 
field due to the radial purely pulsations. In the less stratified sim- 
ulations, such as G8H8, the convective plumes are large enough 
and present an important horizontal filling-factor that leads to a 
velocity field incompatible with the acoustic mode structure. 

Such screening effect, or such spatial distribution of the con- 
vective plumes , have been already found to be responsible of 
the modulation and quenching of gravity modes in DNS of the 
penetrative convection (Dintrans et al. 2005). In our DNS, where 
oscillations are sustained by a physical process acting continu- 
ously, the amplitude of pulsations (and their nonlinear satura- 
tion) may thus be ruled by the extent of the screening of large 
vortices onto the purely radial velocity field. 

6. Conclusion 

In this paper, we have investigated the convection-pulsation cou- 
pling thanks to nonlinear two-dimensional direct numerical sim- 
ulations (DNS). Despite their intrinsic limitations (weak pres- 
sure contrasts across the computational domain or the need of 
large viscosities leading to unrealistic values of Rayleigh's num- 
ber, see e.g. Brandenburg et al. 2000), DNS are without a doubt 
an important way to address this interaction as they fully account 
for the nonlinearities. 

This study follows our former work on the /c-mechanism in 
purely radiative simulations, where we studied both the insta- 
bility and its nonlinear saturation (Gastine & Dintrans 2008a,b, 
Papers I and II, respectively). We have then again modelled 
the ionisation region responsible for the oscillations of classical 
Cepheids by a hollow in the radiative conductivity profile and 
local simulations of a perfect gas layer centered around this ion- 
isation bump have been performed. The initial physical setup has 
been however slightly modified to get a convective zone super- 
imposed with the hollow in radiative conductivity. Furthermore, 
some important numerical improvements have been developed 
(mostly the parallelisation of the ADI solver) as the taking into 
account of strong convective motions coupled with acoustic os- 
cillations requires larger spatial resolutions than the ones used in 
the purely radiative case. 

Using this parallel solver, we have performed 2-D DNS 
of the convection-pulsation coupling in which the oscillations 
are sustained by a continuous physical process based on the k- 
mechanism. Convective motions that develop in these simula- 
tions lead to various impacts onto the acoustic modes: 

- In a first set of DNS, the instability behaves in a similar way 
than observed in purely radiative simulations: there is a lin- 
ear growth and a nonlinear saturation of the acoustic oscil- 
lations. The bulk of kinetic energy is embedded in acoustic 



motions, while convection participates up to 20% to the en- 
ergy budget. However, on the contrary to what was observed 
in Paper II, the nonlinear saturation does not involve any res- 
onance between modes. The frequential analysis emphasises 
numerous harmonics of the fundamental acoustic mode that 
are probably indicative of a mono-mode saturation. In this 
case, the convection does not affect much the oscillations. 
- In a second category of DNS, the convective motions have 
a stronger influence onto the acoustic oscillations. There is 
no clear nonlinear saturation of the /c-instability while the 
amplitudes of acoustic modes remain very weak. Moreover, 
the kinetic energy is almost entirely due to the downward- 
directed convective eddies. The vanishing of harmonics in 
the frequential analysis is a further evidence of the weak- 
ness of the mode amplitudes. In this case, convective plumes 
are quenching the oscillations and this physical phenomenon 
looks like to what is expected to occur in the coldest 
Cepheids close to the red edge of the instability strip. In these 
stars, the large surface convective zone is indeed suspected 
to stabilise the global radial oscillations. 

The various behaviours obtained in these local simulations 
of the /c-mechanism with convection have been further studied 
in the last part of this paper, where we have given some hints to 
explain the physical conditions that lead to the mode quenching. 
Both the influence of the amount of heat carried by convection 
and the density contrast across the layer have been investigated. 
We have first shown that even in the inefficient regime where 
convection only carries 10 - 15% of the total heat flux, the con- 
vective plumes may cancel the radial oscillations, indicating that 
the strength of convection is not the only parameter that explains 
the mode quenching. The density contrast seems on the contrary 
to play an important role on the dynamics of acoustic modes. In 
fact, weaker stratifications (that is, higher values of the density 
scale Hp) correspond to bigger vortices. It means that the scale of 
the more powerful structure in the flow is larger when the strat- 
ification is weaker. In the less stratified simulations, where the 
acoustic modes appear to be strongly quenched by convection, 
the size of convective plumes is indeed larger and that leads to 
an important horizontal filling-factor. In our DNS, the amplitude 
of radial pulsations (and their nonlinear saturation) may thus be 
governed by the screening effect of large convective vortices. 

The 2-D DNS of the /c-mechanism with convection presented 
in this paper are a first step in our study of the convection- 
pulsation coupling occurring in the coldest Cepheid stars with an 
outer convective zone. They are based on a simplified physical 
model that keeps the main physics responsible for the excitation 
of acoustic modes by the /c-mechanism. The opacity bump is for 
instance modelled by a hollow in the radiative conductivity pro- 
file whom the shape is adjusted to get unstable acoustic modes 
superimposed with convection. One of the main prospects of this 
pioneering work is thus to improve this model toward more re- 
alistic stellar conditions (e.g. inputs from tabulated opacities or 
coding of a free- surface boundary condition for the velocity). 
Such developments are however challenging as they will require 
huge numerical improvements in the Pencil Code. 

Another interesting prospect is the comparison between our 
2-D DNS and the prescriptions of different time-dependent con- 
vection (TDC) models (e.g. Stellingwerf 1982; KuhfuB 1986; 
Yecko et al. 1998). In fact, even if our simplified model of the k- 
mechanism is far from the realistic stellar parameters, it could be 
interesting to compare in our simulations the temporal evolution 
of the convective and kinetic fluxes modulated by acoustic oscil- 
lations with their TDC counterparts. We will check the validity 
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of such approaches and give some constraints on the multiple 
dimensionless a coefficients involved in these 1-D models. 
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